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from weak coupling to unitarity: Study of vortices 
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We introduce a nonlinear Schrodinger equation to describe the dynamics of a superfluid Bose gas 
in the crossover from the weak-coupUng regime, where an^^^ <^ 1 with a the inter-atomic s-wave 
scattering length and n the bosonic density, to the unitarity limit, where a — > +oo. We call this 
equation the unitarity Schrodinger equation (USE). The zero- temperature bulk equation of state of 
this USE is parametrized by the Lee- Yang-Huang low-density expansion and Jastrow calculations 
at unitarity. With the help of the USE we study the profiles of quantized vortices and vortex-core 
radius in a uniform Bose gas. We also consider quantized vortices in a Bose gas under cylindrically- 
symmetric harmonic confinement and study their profile and chemical potential using the USE and 
compare the results with those obtained from the Gross-Pitaevskii-type equations valid in the weak- 
coupling limit. Finally, the USE is applied to calculate the breathing modes of the confined Bose 
gas as a function of the scattering length. 
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I. INTRODUCTION 

For the theoretical investigation of the Bose-Einstein 
condensate (BEC) of an ultra-cold bosonic gas the main 
tool is the mean-field Gross-Pitaevskii equation (GPE) 
[l|, , that is reliable for small values of the gas parame- 
ter an^/'^, where a is the s-wave scattering length of the 
inter-atomic potential and n is the bosonic density. By 
manipulating a background magnetic field near a Fes- 
hbach resonance the s-wavc scattering length a can be 
modified and can reach very large values corresponding 
to a strong atomic interaction [3|. To take into account 
the effect of a large scattering length, a modified GPE 
(MGPE) has been introduced by Fabrocini and Polls [1] 
by using the first two terms in the Lee- Yang-Huang ex- 
pansion Q of the energy of a uniform Bose gas. which 
include terms of the order of Vna^. Within the density 
functional approach of Fabrocini and Polls Q , the lead- 
ing term of expansion gives the GPE while the next term 
is responsible for corrections due to moderate atomic in- 
teraction. These two leading terms have been used to 
calculate beyond-mean-field corrections to properties of 
trapped BECs Q . The contribution terms of the order 
of na^ has also been discussed in the literature 0- 

In this paper we generalize the MGPE (3] by consid- 
ering also the behavior of the bosonic system in the uni- 
tarity limit, where a +oo. Jastrow calculations [1, Q 
suggest that in the unitarity limit the zero-temperature 
bulk chemical potential ^ of the Bose system is given 
by /i = ^^n^/^fi^ /m, where ^ is a constant and m is the 
mass of a single atom. (Bulk chemical potential is essen- 
tially the nonlinear term that appear in the mean-field 
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equation and is related to the energy per particle of the 
system.) This result is independent of the atomic scat- 
tering length. In the present work the equation of state 
of the bulk system is parametrized with a Padc approxi- 
mant by using the first two terms of the Lee- Yang-Huang 
expansion and the Jastrow calculations at unitarity 
@, Q. In this way we obtain a time-dependent highly- 
nonlinear Schrodinger equation, that we call unitarity 
Schrodinger equation (USE). The USE in general form is 
time dependent and can be used to study non-stationary 
dynamics, whereas its time-independent form is appro- 
priate to study stationary states. The USE gives the 
hydrodynamic equations of bosonic superfiuids at zero 
temperature, and enables one to study collective dynam- 
ical properties of the system in the full crossover from 
weak-coupling to unitarity. 

As an application of the USE, here we study the struc- 
ture of quantized vortices in both uniform Bose gas and 
Bose gas under axially-symmetric harmonic confinement 
and compare and contrast the results with those ob- 
tained with the MGPE 3 and GPE. Quantized vortices 
in superfiuids are a manifestation of quantum mechanics 
at the macroscopic level [lo| . Recently quantized vor- 
tices have been observed in rotating ultra-cold atomic 
BECs [ll| and also in atomic Fermi gases in the BCS- 
BEC crossover near a Feshbach resonance [l^l- Using 
the GPE, quantized vortices in a harmonic trap have 
been analyzed by Dalfovo and Stringari [l^. More re- 
cently vortices have been investigated with GPE in var- 
ious problems; for instance, collective modes of a vor- 
tex 14 1, vortex under toroidal confinement (isl. [l6j. vor- 
tex lattice [Ol, vortex with attractive scattering length 
[isl . [l9t , collision dynamics of vortices [2^ , collapse of a 
vortex state [2lj , free expansion of vortices [l^l , and vor- 
tex in Bose-Fermi mixtures [2^. Nilsen et al. [l^l used 
the MGPE [1] to investigate the structure of vortices, and 
found a good agreement between the MGPE and varia- 
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tional Monte Carlo results. Using the USE, we extend 
the study of Nilsen et al. [l^l and find that the radius of 
the vortex core decreases with the increase of the scat- 
tering length. At the unitarity limit it reaches a critical 
minimal radius. The properties of this minimal radius 
depend on the trapping geometry: in the case of a vortex 
in a uniform Bose gas it is a decreasing function of the 
uniform density at large distances; in the case of a vor- 
tex in the harmonic trap it is a decreasing function of the 
total number of atoms. Using the USE, we also calculate 
the radial and axial frequencies of collective oscillations 
from weak-coupling to unitarity in a cigar-shaped trap. 

In Sec. [IT] we introduce the present model and relate 
it to supcrfluid hydrodynamics. Formulation for quan- 
tized vortices in the present model is considered in Sec. 
mil Section IIVI is devoted to vortex in a uniform Bose 
gas, where we numerically study the vortex profile and 
vortex-core radius for different scattering lengths. In 
this case the USE is solved by the fourth-order Runge- 
Kutta method [H, [2^ . The vortex-core radius decreases 
with increasing scattering length and saturates to a con- 
stant value in the unitarity limit. In Sec. |V] we con- 
sider a vortex in a Bose gas under axially-symmetric har- 
monic pancake-shaped confinement by solving the USE 
by imaginary time propa gat ion using the semi-implicit 
Crank- Nicholson rule p6l. l27l. [28t . In this case we study 
the profiles of vortices and corresponding chemical poten- 
tials and compare the results with those obtained from 
the MGPE g and CPE. The interesting feature of the 
results of the USE is that they saturate in the unitarity 
limit as a — > 00, whereas the results of MGPE and CPE 
do not saturate in this limit. The frequencies of collective 
breathing oscillatons in a cigar-shaped trap are consid- 
ered in Sec. IVII Finally, in Sec. IVIII we present the 
concluding remarks. 



II. SUPERFLUID HYDRODYNAMICS AND 
NONLINEAR SCHRODINGER EQUATION 

The zero-temperature collective properties of a dilute 
bosonic superfluid under the external potential f/(r) can 
be described by the quantum hydrodynamic equations 
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where n{r,t) is the local density and v(r, i) is the local 
velocity of the bosonic system. In these zero-temperature 
hydrodynamical equations statistics enters the equation 
of state through the bulk chemical potential n{n, a) and 
the quantum-pressure term — [7i^/(2m-y/ri)]V^-\/rI, which 
is absent in the classical hydrodynamic equations [3l|. 
This hydrodynamic regime is achieved in the limit of very 
large number of atoms N. At zero temperature, the bulk 



chemical potential fi{n^ a) of the system is a function of 
the density n and of the inter-atomic scattering length 
a. For a superfluid the velocity field is irrotational, i.e. 
V A V = , and the circulation is quantized, i.e. 



V • dr = 27rL — 
m 
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where L is an integer (angular momentum) quantum 
number p^ . 

We can introduce the complex order parameter 0, [sj 
^'(r,i) = n(r,t)i/2e''S('-,t) such that the phase S{r,t) of 
the order parameter fixes the superfiuid velocity field 
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so that v(r,t) — {h/m)\/ S(r,t). In this way we can map 
Eqs. (P) and ^ into the following time-dependent highly 
nonlinear Schrodinger equation 
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In general, one can use the hydrodynamic equations ([T]) 
and ([2|), or equivalently Eq. to study the global 

properties of the superfiuid, like the stationary density 
profile, the free expansion and the collective oscillations. 

For a Bose gas the following two leading terms of the 
low-density expansion of the bulk chemical potential can 
be obtained [2J] from the expression for energy per par- 
ticle as obtained by Lee, Yang and Huang 
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where v}/^a is the dimensionless gas parameter js^ . 
Note that in this expansion the scattering length a must 
be positive (a > 0) corresponding to a repulsive interac- 
tion. Higher order correction terms to the bulk chemical 
potential have also been considered in the literature 0]. 
The lowest order term of expansion ^ was derived by 
Lenz [s^. Considering only this term, Eq. ([5]) becomes 
the familiar CPE [J, Q 
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By taking into account also the second term of the ex- 
pansion ([6]), Eq. (O becomes the so-called modified 
Gross-Pitaevskii equation (MGPE) introduced by Fab- 
rocini and Polls H: 
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In the unitarity limit, where a — > -l-oo, for dimensional 
reasons [1, Q the bulk chemical potential must be of the 
form i 



u(n, a) = £ — n 
m 



2/3 
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where ^ is a universal coefficient. Thus, in the unitarity 
limit the bulk chemical potential is proportional to that 
of a non-interacting Fermi gas [s^j- Recent numerical 
calculations based on Jastrow variational wave functions 
give ^ = 22.22 and we shall consider this value of ^ in 
our calculation. 

In the full crossover from the small-gas-parameter 
regime to the largc-gas-parameter regime, we suggest the 
following expression as the bulk chemical potential of the 
Bose superfluid 



(10) 



where f{x) is an unknown dimcnsionlcss universal func- 
tion of the gas parameter x = n^/^a. A general Pade 
approximant for the function f{x) consistent with the 
expansion of Lee, Huang, and Yang ^ Q and Eq. fTU]) 
should have the following form 



f{x) = 47r. 
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where a, (3, and 7 are yet undetermined parameters. 
Without further information about /i(n, a) we cannot de- 
termine all these parameters consistently. In this paper 
we consider the minimal form of this function consistent 
with Eqs. ^ and pI7|) obtained by setting 7 = and a = 
32/(30F) and (3 = 47ra/^, with i ^ 22.22. The function 
/(x) of Eq. (ini) is such that f{x) = 4,tt[x+32x^^^ /{3^)] 
for a; ^ 1 and the present model reduces to the MGPE 
([5]) [I]. In the opposite extreme x ^ 1, f{x) = and 
the present model reduces to the unitarity limit ([9]). We 
call Eq. ^ equipped with Eqs. pUj) and (dJ), e.g.. 
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by the name USE, i.e. unitarity Schrodinger equa- 
tion. As, by construction, the USE has the proper weak 
[23 | and strong Q coupling limits, it is appropriate for 
the study of weak-to-strong coupling crossover. This is 
the model equation we use in the following sections to 
study quantized vortices in a Bose superfluid from weak 
to strong coupling. In this Sec. 5* is normalized by 
J\^\^dr = N. 

It is important to stress that the present approach, 
based on the quantum hydrodynamics, could be applied 
to both superfluid bosons and fermions [l^, [IHl- This 
hydrodynamic approach is a time-dependent local den- 
sity approximation with gradient corrections. In the last 
few years we have succesfully applied it to investigate 
the collective properties of different dilute systems, like 
Bose-Fermi mixtures 13611 , the superfluid Fermi gas in the 
BCS-BEC crossover [331, H-* Lieb-Liniger liquid 

[3^. Here we have obtained the USE valid from weak- 
coupling to unitarity based on the general quantum hy- 
drodynamical scheme. 



III. QUANTIZED VORTICES 

The structure of vortices in superfluid ''He at zero- 
temperature was investigated many years ago by Chester, 
Metz and Reatto [s^ by using a many-body variational 
wave function, and more recently by Dalfovo (40| by us- 
ing the Orsay-Trento density- functional [4l|. The two 
approaches, which give very similar results, are based on 
the assumption that the superfluid velocity of the quan- 
tized vortex rotating around the cylindric z axis is given 

by 
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where L is the quantum number of circulation, p is 
the cylindric radial coordinate and is the unit az- 
imuthal vector, with </) the azimuthal angle. In 1997 
Sadd, Chester and Reatto [i^] suggested that the veloc- 
ity of superfluid "^He is not truly singular at the vortex 
line and the vorticity is distributed over a flnite region. 
Nevertheless, the deviations from Eq. in ^He seem 
to be very small [i^ . 

We recall that there are remarkable differences between 
superfluid '^He and the dilute superfluid we are consid- 
ering here. In *He the effective radius Rq of the inter- 
atomic potential is of the order of the average distance 
n~^/^ between atoms, i.e. Rqv}/^ ~ 1. Instead, in dilute 
ultra-cold gases the effective radius i?o is always much 
smaller than the average distance n~^/^ between atoms, 
i.e Rqv}/^ <^ 1. For a dilute gas the coupling regime de- 
pends on the scattering length a, namely on the gas pa- 
rameter an^^^: in the weak-coupling regime an^^^ ^ 1, 
while in the strong-coupling regime an^^^ ^ 1. Another 
remarkable difference between liquid helium and quan- 
tum gases of alkali-metal atoms is the width of the vor- 
tex core. The core of a quantized vortex is only few 
angstroms in superfluid '*He while it is of the order of 
sub- micron in a dilute atomic gas [43] . 

In our USE we get Eq. (jTS]) by setting *(r,<) = 
Tpip, z) exp[i(L0— ^^)] ' where /ip is the chemical poten- 
tial of the inhomogeneous superfluid, fixed by the nor- 
malization. In this way, Eq. Q with Eq. (jTU]) becomes 
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+ U{p, z)+ p.{n,a) ^ = fiQ 
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i^(p^) is the Laplacian operator in the 



radial direction and 



" is the centrifugal term which 
determines the size of the vortex core, that is of the order 



of the heahng length Ih = fiL / ^2mp,Q [4^. This expres- 
sion is obtained by equating the centrifugal term to the 
chemical potential /ip. 

We introduce scaled variables by using the character- 
istic length Ic of the system. In particular, we make the 
following transformations: p = p/lc, z = z/lc, a = a/l^. 



U = U{mll)lh\ po = Po{mll)lh\ V 
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FIG. 1: (Color online) Numerical results for scaled density of 
vortices with L = 1 and 2 in a uniform Bose gas ob- 

tained by solving the USE for scaled scattering lengths a/lc = 
(a) 0.01, (b) 0.03, (c) 0.1, and (d) 10 compared with the an- 
alytic asymptotic results (|20|l or ((21}. 



IV. VORTEX IN A UNIFORM BOSE GAS 

Let us first consider the case U{p, z) = 0. We suppose 
that the Bose gas is asymptotically uniform, i.e. 



1 
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for p,z — > oo. This condition fixes the characteristic 

— 1/3 

length Ic, that must be Ic = noo with the uniform 
density at infinity. In this case the chemical potential is 
fixed by the asymptotic condition, from which one finds 



Mo = ./(a) 



(17) 



In addition, we can choose a wave function which depends 
only on p, and Eq. p5|) becomes 



1 Id 
29^ ~ 2^9^ ^ 2^ 



= . 

(18) 

This second-order ordinary differential equation must be 
solved numerically. Clearly with i ^ the vortex func- 
tion i/)(yo) has a core around p = 0. The asymptotic 
behavior olip^p) for p ^ oo is obtained by neglecting the 
spatial derivatives in Eq. (jlSp . In this way we obtain the 
algebraic equation 



m^/'frnf^o^) = /(«) 



2p2 



(19) 



In the weak-coupling regime, a <^ 1 where f{x) = Attx, 
Eq. (HH) gives 



fim Hp) = ( 1 - -4^ 



1/2 
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Instead, in the strong-coupling regime, a 3> 1 where 
fix) = ^, Eq. m gives 



lim V(p) = ( 1 - 777^ 
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Although the limit (PO)) depends on the scaled scattering 
length a, the limit ((2T|) is independent of a. The scaled 
healing length Ih = Ih/lc, which estimates the size of the 
vortex core [3l , is given by 



way Eq. (fT4|) becomes 
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Now we consider Eq. (fT5|) in the limit p ^ 0. Because 
of the regularity of the wave function and the function 
f{x) in this limit the last two terms [the terms involv- 
ing the function f{x)] in this equation remains finite and 
can be neglected in comparison to the angular momen- 
tum term L'^ / {2p^). Consequently, we have the following 
condition at small p 



where '0(Pi 2) is assumed to be real. 
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FIG. 2: (Color online). The numerical values of vortex core 
radius for L = 1 and 2 vs. scaled scattering length calculated 
by the USE compared with the analytic results of healing 
length given by Eq. . 



TABLE I: The parameters for L = 1 used in the numerical 
solution of Eq. (|18|l as plotted in Fig. \T\ in this case 1/5(0) = 0. 
For L = 2, ^(0) = ij'{0) = 0. 
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FIG. 3: (Color online). The numerical values of effective 
nonlinearity N"^^^ f{N^/'-^a) of the USE, 4nNa of the GPE, 
and 47rAfa(l + aa^^^VN) of the MGPE for i/> = 1 for (a) 
N = 1000 vs. a/lc and (b) a/lc = 0.05 vs. iV. 



where f{p) is a smooth function. 

We next solve Eq. (fT8|) for L = 1 and 2 using the 
classic fourth-order Runge-Kutta method [1^, [2^ . This 
method produces very accurate convergence. We employ 
a space step of A = 0.0001 and integrate up to a scaled 
distance of prnax = 20. The integration is started with 
the initial boundary condition (|23p with a trial '0(0) and 
■0'(O). For L = 1 this condition is taken as ■0(0) = and 
ip'{0) = constant and the integration started at p = 0. 
For L = 2 both "0(0) = and '(/''(O) = and the numerical 
integration cannot be started at p = 0, as then ijj as well 
as 0' at subsequent sites become zero. For L=2, the inte- 
gration is started at p = A with the boundary condition 
■0(A) = and ?/''(A) equal to a small constant. The inte- 
gration is propagated to p = prnax where the asymptotic 
conditions ((20)l or ((2T|) remain valid. If after integration 
with a trial guess, the boundary conditions or (PT]) 
cannot be satisfied, the method is implemented with a 
new trial guess. The process is continued until a solution 
satisfying the proper boundary conditions at small and 
large p is obtained. We obtain the solution for different 
values of a. 

The numerical results for the scaled density [= '0^(p)] 
for different scaled scattering length a = a/lc is plotted 
in Fig. [1] and compared with the asymptotic forms ([20]) 
or (j2ip for L = 1 and 2. The corresponding parameters 
■0'(O) and /(a) used in the numerical solution arc given in 
TablclH The quantity /(a) shown in Tablc|T]is interesting 



as it can be compared with the chemical potential of Eq. 
(fT8)) . Hence the numerical values of /(a) for different a 
should give an idea of the variation of chemical potential 
with coupling. The quantity /(a) is also related to the 
scaled healing length Ih given by Eq. ((22)) as shown in 
Table H 

It is interesting to calculate numerically the vortex core 
radius defined as the value of p for which ip'^{p) = 0.5. In 
Fig. [5] we plot the numerical values of vortex core radius 
versus scaled scattering length and compare with the cor- 
responding analytical results of healing length given by 
Eq. (P^ . The figure shows that the vortex-core radius 
has the same trend of the analytical healing length [i^ ]. 
From Figs. [T] and [2] we find that the vortex-corc radius 
decreases with increasing coupling but increases with an- 
gular momentum L. However, for very large coupling it 
saturates. This saturation of the vortex core is properly 
given by the USE. 



V. VORTEX IN A BOSE GAS IN A HARMONIC 
TRAP 

Now we consider the bosonic fluid in an axially- 
symmetric harmonic trapping potential, i.e. 

U^^mcvlip' + X^z^) (24) 
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where A = Wz/ujp is the anisotropy parameter, with ujp 
the transverse frequency and ojz the axial frequency. We 
choose as characteristic length of the system the trans- 
verse harmonic length, i.e. Ic = Ti / {muj p) . By us- 
ing the scaled variables the confining potential reads 
U — -I- A^z^) . With this confining potential the 

mean-field equation can be written explicitly as 



1 + aa3/2iVi/2^ 



(25) 



l-f/3a5/2iV5/6V;5/3_ 

Now the normalization condition is given by 

2tt P dp dz ■ijj {p,z) = 2tt p dpricip) = 1(26) 
Jo J -oo Jo 

Here we defined ndp) as the column density (as in [23 |) 
appropriate for the study of radial distribution of mat- 
ter in the vortex. Equations ((25)) and (|26|) constitute the 
USE for an axially trapped Bose gas valid from weak- 
coupling (a — > ) to the unitarity limit (a —^ oo) and 
(in spite of a slightly complicated numerical structure) 
it is no more difficult to solve than the usual GP equa- 
tion valid in the weak-coupling limit. In the extreme 
weak-coupling limit, the non- linear term of Eq. (|25p be- 
comes the usual GP term AirNaip'^ . In the unitarity limit 
this term becomes independent of scattering length and 
equals ^7V2/3?/;4/3, ^ = 22.22. We recall that if we use the 
GPE dll), as well as the MGPE (O, in the unitarity hmit, 
it will lead to an inappropriate dependence of the non- 
linearity on scattering length, as well as on the number 
of atoms N. 

It is appropriate to study the behavior of the nonlin- 
earity of the USE that is iV^/a^Vs j(jvi/3^2/3^^^ 

To study the dependence of this nonlinearity on N and 
a we set ?/; = 1 in this expression. The resultant nonlin- 
earities are plotted in Fig. [3] (a) and (b) as a function 
of a and N for constant N = 1000 and a/lc = 0.05, re- 
spectively. The interesting feature of the nonlinearity of 
this equation is exhibited in Fig. [3] (a), where we find 
that, for a fixed TV, the nonlinearity of this equation sat- 
urates at large scattering lengths a. In the GP model 
the nonlinearity increases linearly for all a, whereas for 
the MGP model it increases indefinitely, but with a more 
complicated dependence on a. From Fig. [3] (b) we find 
that for a fixed scattering length the nonlinearity of Eq. 
(PS)) as well as the GP equation increases with N. For 
the GP equation it increases linearly with N, whereas for 
Eq. ([25]) it increases as A^^/s j-qj. ^^j-gc 

To further study the USE (|^5|) , we solve it numerically. 
For numerical convenience we transform this equation to 
its time-dependent form by replacing the term pgip by 
its time-dependent counterpart idip/dt appropriate for 
non-stationary states. However, in this paper we limit 
ourselves to a study of stationary states. Nevertheless, 
the introduction of time-dependence allows for a solution 
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FIG. 4; (Color online). The numerically calculated scaled 
chemical potential p,o (in units of oscillator energy Tiujp) from 
GPE, MGPE and USE vs. (a) scaled scattering length a/U for 
N = 10000, A^ = 8, and vs. (b) N for a/U = 0.15155, A^ = 8. 



of Eq. (j25p by imaginary time propagation method af- 
ter discretizing it with the semi-implicit Crank-Nicholson 
rule [2^, [13, HI] . In the process of discretization we use a 
time step of 0.0005 and a space step of 0.03. We limit our 
numerical study to the L = 1 vortex, as i > 1 vortices 
are unstable and decay to several L = 1 vortices conserv- 
ing the angular momentum. In the numerical simulation 
of a vortex with L = 1 we employ a pizza-shaped con- 
densate with trap anisotropy A^ = 8 as in the theoretical 
study of Nilscn et al. [l^l and the experiment at JILA 
[4^. After the wave function '4}{p) of Eq. ([^5]) is obtained 
by the imaginary time propagation, the chemical poten- 
tial is obtained by multiplying this equation by '0(p) and 
integrating over all space. 

In Fig. |4] (a) we plot the numerical results of scaled 
chemical potential vs. scaled scattering length a/lc as 
obtained from the USE ^ as well as the GPE 
for N = 10000. In Fig. H (b) we plot the same vs. 
N for a/lc = 0.15155. The scaled scattering length 
a/lc = 0.15155 is the value studied by Nilscn et al. and 
is appropriate for Rb atoms in an experimental trap used 
at JILA for a scattering length a = 35a(Rb), where 
a(Rb)= lOOao (ao the Bohr radius) is the experimen- 
tal atomic scattering length of Rb. In Table [IT] we plot 
some of the results for chemical potential /2o for different 
a and N values. The first row in this Table agrees with 
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TABLE II: Scaled chemical potential ^io of GPE, MGPE, and 
USE for different a — a/l^, N, and GP nonlinearity dN. 
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FIG. 5: (Color online) A typical profile of the probability 
density P{p/lc, z/lc) — ip^ {p/lc, z/lc) for a vortex obtained 
with the USE for N = 10000 and a = 1. 
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FIG. 6; (Color online). The numerically calculated column 
density nc{p) = J^_^dzijj'^{p,z) from the USE, GPE, and 

MGPE vs. p/h for (a) = 10000, = 8 and different 
a/lc, and (b) a/lc = 0.15155, = 8 and different N. 



the calculation of Nilscn et al. for the GP and the MGP 
models. 

From Fig. |4] (a) and Table |TT] we find a saturation of 
the present mean-field result for large scattering length 
in the unitarity limit. This is consistent with the satu- 
ration of the nonlinearity of the USE ([^5]) in this limit. 
After this saturation is obtained, any further increase in 
the scattering length at a fixed N does not change the 
chemical potential /io . Figure S] must be compared with 
Fig. [3l For a fixed scattering length, as N is increased 
the nonlinearity of the USE (QH) keeps on increasing. As 
a consequence the chemical potential also increases indef- 
initely as N . In Fig. 2] (b) the USE chemical potential is 
always greater than the GP one. On the other hand, for 
a fixed N , as a is increased the nonlinearity of the USE 
saturates above a certain value of a. This has a con- 
sequence in the results: With the increase of scattering 
length, as the system tends to the unitarity limit, there 
are a saturation of the present chemical potential and a 
crossover, beyond which the USE chemical potential be- 
comes smaller than the GP one. From Figs. [Hwe find 
that for moderate to small a values the results for chem- 
ical potential of the MGPE ^ 0] and USE (USD remain 
very close to each other. From Table HIl we find that for 
N = 10^, and a = 0.01 (corresponding to a large GP 



nonlinearity of Na = 10000) the chemical potentials of 
the MGPE and the USE differ by less than 0.2% whereas 
the chemical potentials of these two models differ from 
that of the the GP model by about 3%. (We recall that 
a = 0.01 is the typical experimental value of this quan- 
tity in the experiment at JILA 0, [i^.) The chemical 
potentials of MGPE and USE start to differ for large a 
values (larger than the experimental value but attainable 
by the Feshbach resonance technique Q): the chemical 
potential of the USE model exhibits saturation, whereas 
the chemical potential of the MGP model increases very 
rapidly with increasing a values which has been made ex- 
plicit in the data in the last two rows of Table |TT1 where 
the results of the MGPE 0] and USE show the main 
differences. 

In Fig. [5] we demonstrate a typical profile of the prob- 
ability density of a vortex with L = 1 as obtained from 
the USE for N = 10000 and a = 1, corresponding to a 
large GP nonlinearity of Na = 10000. The density is 
manifestly zero for p = for the vortex state. However, 
both for theoretical and experimental analysis it is ap- 
propriate and more convenient to calculate the column 
density nc{p) as defined in Eq. ([26]) to study the ra- 
dial distribution of matter in a vortex for the USE, GPE 
and MGPE. The results for column density are plotted 
in Figs. [Hfor (a) N = 10000 for different l/lc and for 
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FIG. 7: (Color online). The numerically calculated scaled rms 
sizes (p/lc) and {z/lc) of the condensate from USE, GPE, and 
MGPE vs. scaled scattering length a = a/lc 



(b) a/lc = 0.15155 for different N. The densityprofile 
of Fig. [S] (b) for N = 500 is also reported in [2J| and 
the two agree with each other for the GP equation. An 
interesting feature of the plots of Figs. [S] is that with 
an increase of nonlincarity (cither via an increase of a/lc 
for a fixed A^, or via an increase of N for a fixed a/lc), 
the column density extends up to a larger radius. Of the 
column densities, that obtained with the USE extends to 
a larger radius than the GP column density in general. 
This trend is reversed in the unitarity limit of large scat- 
tering length as can be seen from the plot of Fig. [S] (a) 
for a/lc = 3. For small values of a/lc the MGPE results 
are close to the GPE results. For medium values of a/lc 
the MGPE results are close to the USE results. However, 
for large values of a/lc the MGPE wave functions extend 
to very large distances, compared to the other models, 
due to an unphysically large nonlincarity in this model 
(MGPE). 

A manifestation of the saturation of nonlincarity of 
the USE is explicit in the rms (root mean square) ra- 
dial {{p/lc)) and axial {{z/lc)) sizes of the condensate as 
the scattering length a is increased. The results for the 
rms sizes calculated with the GPE keeps on increasing 
indefinitely as the scattering length is increased. The 
rms sizes calculated with the MGPE increases even more 
rapidly than the results of the GPE. However, the rms 
sizes calculated with the USE saturates after a certain 
value of a, beyond which they remain constant. This 
is illustrated in Fig. [7] where we plot the rms sizes vs. 
log(a/;c) as calculated by the USE, GPE, and MGPE. As 
we are considering a pizza-shaped condensate the radial 
size {p/lc) is larger than the axial size z/lc- The theoret- 
ical results for the rms sizes exhibited in Fig. [7] can be 
verified experimentally with present technology and this 
will provide a test for the USE proposed in this paper. 
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FIG. 8: (Color online) Square of the reduced frequencies of 
collective radial and axial oscillations (a) Up and (b) in a 
cigar-shaped trap vs. atan(a/Zc). 

VI. COLLECTIVE OSCILLATIONS IN 
HARMONIC CONFINEMENT 

We consider the effect of confinement due to an exter- 
nal anisotropic harmonic potential (|24p on frequencies of 
collective oscillation of the system from weak-coupling to 
unitarity using the USE. It has been shown by Cozzini 
and Stringari [46| that assuming a power-law dependence 
fi = AnF for the chemical potential (polytropic equation 
of state [13]), from Eqs. ([1]) and ^ without the quan- 
tum pressure term, one finds analytic expressions for the 
collective breathing frequencies of the superfluid. In par- 
ticular, for the very elongated cigar-shaped traps, the 
collective radial breathing mode frequency f2p is given 
by 

ftp = V2(r + 1) LUp , (27) 
while the collective longitudinal breathing mode fl^ is 

Here wc introduce an effective polytropic index P as 
the logarithmic derivative of the bulk chemical potential 
/i, that is 

ndn 2 1 f'{x) , , 

fi on 3 3 j(x) 
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where x = ri^/^a is the gas parameter. This formula for 
the local polytropic equation is useful to have a simple 
analytical prediction of the collective frequencies. In the 
weak-coupling regime (x <^ 1) one finds xf'{x)/f{x) ~ 1 
and r = 1, consequently, the reduced frequency square 
ly^ = (Tlp/iVpy = 2(r + 1)) = 4 and = f = 

3 — l/(r-|-l) = 5/2; while in the unitarity regime {x 3> 1) 
it holds xf'{x)/f{x) = and F = 2/3 and consequently, 
i^p = 10/3 and z/^ = 12/5 . The analytical predictions 
of Eqs. dSZl) and ^ with Eq. ^ are shown in Fig. 
ISl where we plot i^^ and v'^ vs. atan(a/^c), so that the 
entire region oo > a > is mapped into the finite interval 
tt/2 > atan{a/lc) > 0. 

VII. CONCLUSION 

In this paper we have proposed a new nonlinear 
Schrodigcr equation — Eq. (fT2|) or ([25]) — to study the 
properties of a BEC or a general superfluid Bose gas valid 
in both the weak-coupling and strong-coupling (or uni- 
tarity) limit. We call this equation the USE (unitary 
Schrodiger equation). This equation has a complicated 
nonlinearity structure in its dependence on scattering 
length a and number of atoms N. In the extreme weak- 
coupling limit the USE reduces to the usual mean- 
field GP equation O 0; for medium coupling it becomes 
the modified GPE (MGPE) © introduced by Fabrocini 
and Polls Q , which is a generalization of the GP equation 
valid for medium coupling incorporating the correction to 
the bulk chemical potential of a Bose gas as introduced 
by Lee, Yang, and Huang Q. However, for very strong 
coupling, both the GP and the MGP equations break 
down and the bulk chemical potential attains a satura- 
tion Q . Considering this bulk chemical potential Cowell 
et al. @ derived a nonlinear equation for a superfluid 
Bose gas valid in the strong coupling limit. The USE 
reduces to the equation by Cowell et al. in the strong 
coupling limit and thus have the correct form in both the 
weak and strong-coupling regimes. Hence this equation 
should be useful to study the crossover of a superfluid 
Bose gas from the weak to strong-coupling limits. In the 



time-independent form the USE is useful to study the 
stationary properties of a BEC. The full time-dependent 
USE can be used to study non-equilibrium properties of 
a dilute BEC, such as, dynamical oscillations [2^, col- 
lapse [11] , free expansion after release from the trap [l^l , 
soliton formation and soliton dynamics [13, etc. 

In this paper we applied the USE to study vortices in a 
superfluid Bose gas. First, we study vortices in a uniform 
Bose gas as the scattering length is varied from weak to 
strong coupling values. The vortex core radius decreases 
with increasing scattering length eventually attaining a 
saturation value. The vortex core radius is comparable 
to the healing length in this case. 

Next, we study vortices in a BEC confined by an 
axially-symmetric harmonic trap by solving the USE nu- 
merically. It is found that the effective nonlinearity of 
the USE saturates with the increase of scattering length 
a in the strong-coupling limit for a fixed number of par- 
ticles N. In this limit the chemical potential as well as 
the rms sizes of the BCS saturate attain constant values. 
However, the rms sizes obtained from the USE and the 
MGPE grows indefinitely as the scattering length is in- 
creased towards the unitarity limit. In the weak-coupling 
limit the results of the USE are compatible with those ob- 
tained from the GPE and the MGPE. Finally, we present 
results of axial and radial frequencies of collective oscil- 
lation of a BEC in a cigar-shaped trap using the USE. 

The results of this paper for a BEC confined in a 
axially-symmetric harmonic trap can be tested experi- 
mentally with present technology specially in the uni- 
tarity limit near a Fcshbach resonance and this will 
provide a stringent test for the proposed USE. 
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